Faster but Not Sweeter: A Model of Escherichia coli Re-level Lipopolysaccharide for Martini 3 and a Martini 2 Version with Accelerated Kinetics

Lipopolysaccharide (LPS) is a complex glycolipid molecule that is the main lipidic component of the outer leaflet of the outer membrane of Gram-negative bacteria. It has very limited lateral motion compared to phospholipids, which are more ubiquitous in biological membranes, including in the inner leaflet of the outer membrane of Gram-negative bacteria. The slow-moving nature of LPS can present a hurdle for molecular dynamics simulations, given that the (pragmatically) accessible timescales to simulations are currently limited to microseconds, during which LPS displays some conformational dynamics but hardly any lateral diffusion. Thus, it is not feasible to observe phenomena such as insertion of molecules, including antibiotics/antimicrobials, directly into the outer membrane from the extracellular side nor to observe LPS dissociating from proteins via molecular dynamics using currently available models at the atomistic and more coarse-grained levels of granularity. Here, we present a model of deep rough LPS compatible with the Martini 2 coarse-grained force field with scaled down nonbonded interactions to enable faster diffusion. We show that the faster-diffusing LPS model is able to reproduce the salient biophysical properties of the standard models, but due to its faster lateral motion, molecules are able to penetrate deeper into membranes containing the faster model. We show that the fast ReLPS model is able to reproduce experimentally determined patterns of interaction with outer membrane proteins while also allowing for LPS to associate and dissociate with proteins within microsecond timescales. We also complete the Martini 3 LPS toolkit for Escherichia coli by presenting a (standard) model of deep rough LPS for this force field.


Details for setting up the atomistic (CHARMM36m) simulation
NPT conditions were maintained using a Nosé-Hoover thermostat 1,2 , and a Parrinello-Rahman barostat 3 with a compressibility of 4.5 × 10 −5 bar −1 .The integration timestep was set to 2 fs, and bonds to hydrogens were constrained with the LINCS algorithm 4 .Long-range electrostatics were computed using particle-mesh Ewald 5 with a tolerance of 10 −6 , 4th-order spline interpolation, and a maximal mesh size of 0.12 nm; van der Waals interactions were shifted smoothly to zero between 1.0 and 1.2 nm.

Details for setting up coarse-grained (Martini 2 and 3) simulations
All the systems were minimized using the steepest descent algorithm until the forces converged to less than 1000 kJ/mol/nm.With the Martini 2 force field, the systems were equilibrated once at 313 K using a 20 fs timestep, except for the PMB1 simulations where a 10 fs timestep was used.The time of equilibration for simulations of the OMP island and PMB1 were 30 ns and 100 ns respectively.All other Martini 2 systems were equilibrated for 5 ns each.During equilibration, the cut-off algorithm combined with the potential-shift modifier was employed to treat both van der Waals and coulombic interactions.
The van der Waals and coulombic radius was set to 1.1 nm.The pressure was maintained at 1 bar with semi-isotropic pressure coupling using the Berendsen barostat 6 , the compressibility was set to 3.0-4 /bar, the temperature was maintained at 313 K using the V-rescale thermostat 7  For steered molecular dynamics (MD) simulations, this equilibrated snapshot was used as the initial frame and a force of 1000 kJ/mol/nm was applied on the center of mass of DPC in the -Z direction to pull it from water into the bilayer using at the rate of 1 nm/ns.These steered MD simulations were performed in 5 replicas using the standard Martini 2 and Martini 3 ReLPS models.All production simulations including the steered MD simulations used the Parrinello-Rahman barostat 3 along with the velocity rescaling thermostat 7 .

Simulation analyses and visualization
The GROMACS package 8 version 2019.4 or higher was used to setup, run, and post-process all simulations.PyLipID 9 was used to calculate residence times between ReLPS and outer membrane proteins (OMPs) using a lower and upper cutoff of 0.485 and 0.8 nm to define contacts.RMSD, distances, lateral diffusion, density maps, radial distribution functions and contacts were analyzed using GROMACS analysis tools: gmx rms, gmx distance, gmx msd (-lateral z), gmx densmap, gmx rdf, gmx mindist respectively.The lateral diffusion coefficients of the Martini 2 ReLPS models were obtained for independent 1 µs-long segments from a linear fit of the MSD using a time window of 20-100 ns for the fit.Area per lipid (APL) and membrane thickness were analyzed using FATSLiM 10 .Phosphate moieties were used to define the group for leaflet identification.Lipid order parameters were calculated using a modified version of the do-order-multi.pyscript available on the www.cgmartini.nlwebsite (http://www.cgmartini.nl/index.php/downloads/tools/229-do-order).VMD was used for visualization of all systems 11 .NumPy 12 was used during postprocessing and analyses of data.Matplotlib 13 and Xmgrace (https://plasma-humphreygate.weizmann.ac.il/Grace/) were used for plotting.The nearneighbour (NN) index 14 was calculated using the following equation:

NN Index = (unique neighbour number -NNN) / (NLPS -1 -NNN)
where unique neighbour number refers to the number of unique LPS molecules contacted by the reference LPS throughout the trajectory, NNN is the number of surrounding LPS molecules in a single frame of the trajectory (taken as 6), NLPS is the number of LPS molecules in the simulation.

Details of obtaining the fast ReLPS Martini 2 model
We first investigated reduction of the magnitude of the charges (by 100, 50 and 25 %) on all anionic beads of ReLPS.Each of these settings were tested with neutralizing monovalent (Na + ) and divalent (Ca 2+ ) cations.To quantify the impact of these alterations to charges we consider the following metrics: Concomitantly, the PO-PO RDF curves reveal the reverse trend with decreasing charge.As more and more POs come closer to each other, they repel due to electrostatics leading to increased MSD and NNI.However, at the point at which the charge is totally removed from the PO bead, despite their proximity, there is no electrostatic repulsion between PO-PO beads and therefore the average MSD and NNI drop slightly.We provide the non-bonded potential plots for each charge regime to further highlight this point (fig S7A).The same trend is observed with PO-Na + interactions (fig S7B).
The partial density profiles of various components of this charge-modified (Q50%, Na + ) system suggest that the system is stable and comparable with those of the standard Martini 2 system (fig S8A).
Concomitantly, there was a minimum increase in the APL from 1.67 ± 0.23 nm 2 to 1.73 ± 0.3 nm 2 (fig

S8B
).Therefore, based on these data, the model whereby the charges of both the ReLPS anionic beads and neutralizing cations were scaled down by 50% was deemed to provide the best compromise between deviation from the standard Martini 2 model and faster kinetics, and was taken forward to the next step.
We next investigated scaling down the inter-ReLPS LJ interactions (which leads to a shallower welldepth of potential energy).We tested a variety of scaling combinations i.
ReLPS CHARMM36m (all-atom) to Martini 3 (CG) mapping scheme.In (A), the dark line represents the mean value from the three trajectories and the shaded regions indicate the standard deviation.Systems are labelled according to the modifications made to the anionic charges of ReLPS beads and the ions used to neutralize the system.For example, "Q50%, Na + " means that the negative charges on all anionic ReLPS beads were scaled down to 50% and the system was neutralized by monovalent cations in addition to 0.15 M NaCl.The MSD and NNI data are split into two plots each for visual clarity.

S1. Martini
and a time constant of 1 ps.With the Martini 3 force field, systems were equilibrated step-wise using the Berendsen barostat for a total time of 150.6 ns i.e. (1) with 1 fs timestep at 200 K for 50 ps, (2) with 1 fs timestep at 313 K for 50 ps, (3) with 5 fs timestep at 313 K for 500 ps, (4) with 10 fs timestep at 313 K for 50 ns, (5) with 10 fs timestep at 313 K for 100 ns.From step 3 to step 5, the protein backbone and ReLPS phosphate beads were position restrained with 50 kJ/mol/nm harmonic potential.For equilibration and production stages, the reaction-field algorithm was used to treat coulombic interactions and the cutoff with the potentialshift modifier algorithms to treat van der Waals interactions in Martini 3 simulations.For the membrane permeability simulations, the octane/ DPC/ DOPC molecule was position restrained during equilibration.
e. scaling the LJ interactions of only GlcN headgroup beads, only the Kdo core (S beads), and both GlcN and Kdo core beads simultaneously (fig S9A-C).The scaling regimes were evaluated by calculating MSD and NNI of ReLPS.These data revealed that scaling only the Kdo core had little impact on MSD and NNI, whereas greater impact was achieved when modifying either only the GlcN headgroup beads or GlcN and Kdo core beads.Thus, in the interests of limiting chemical modifications as much as possible, only alterations to the GlcN beads were considered hereon.Compared to the standard Martini 2 model, while scaling down the LJ interactions of the GlcN beads by 5 levels (as defined by Martini 2 conventions) led to the greatest lateral diffusion (fig S9B), we took into account the (1) change in the APL i.e. an increase less than 5% (table S3), (2) the partial density profiles showing a good agreement with the standard Martini 2 model (fig S10A), and (3) membrane thickness showing an increase (~1%) less than 0.1 nm (fig S10B).As a result, we considered the model with LJ scaled down by 2 levels.

3 1 S2. Martini 3
ReLPS mapping and associated bead types.Groups of two to four heavy atoms were mapped into coarse-grained (CG) beads.Notably the sugar rings were modelled with an extra virtual site bead (depicted in green) to better characterize their volume and interactions.Here, each highlighted segment of the 2D depiction corresponds to the same-coloured bead on the right CG molecule taken from a simulation snapshot.The four different colours represent: Kdo sugars i.e. first two in core region (yellow), GlcN sugars from lipid A moiety (grey), phosphate substituents (brown), acyl chains (cyan).The labels on the CG representation show the final assigned beadtype.ReLPS bond distributions.Distributions for all the bonds are plotted separately.Blue bins correspond to the Martini 3 model whereas black outlines correspond to the data obtained from pseudo-coarsegrained atomistic simulation (i.e converting the atomistic coordinates into coarse-grain model for each of the frames).S3.Martini 3 ReLPS angle distributions.Distributions for all the angles are plotted separately.Blue bins correspond to the Martini 3 model whereas black outlines correspond to the data obtained from pseudo-coarse-grained atomistic simulation (i.e converting the atomistic coordinates into coarse-grain model for each of the frames).its (A) MSD and (B) NNI.For each plot, the data shown are calculated from three simulation replicas.

S7.
Reasoning for MSD and NNI of ReLPS approaching a peak and then decreasing when their phosphate bead charges are scaled from 100% to 50%, 25% and 0%.Total potential energies (LJ + coulombic potentials) (top row) and RDFs (bottom row) of monovalent, divalent cations, and other LPS phosphates (PO beads) with reference to LPS phosphates, labelled as 'PO-NA', 'PO-CA', 'PO-PO' respectively.The system charge neutralizing cations are (A) the divalent cations (Ca 2+ ), and (B) the monovalent cations (Na + ).S8. Evaluating the influence of only scaling down charges on Martini 2 ReLPS beads and cations on its partial density profiles and area per lipid.(A) Partial density profiles of various components of the 'Q50%, Na + ' system (top half) compared to the (unmodified) reference system 'Q100%, Ca 2+ ' (bottom half).The density profiles are such that the hydrophobic core of the membrane is located at the i.e.Z distance = 0 nm (grey dotted line), and the periplasmic and extracellular regions are located on the left and right hand side of the dotted line respectively.(B) APL versus time averaged across simulation replicas of the various systems.The reference system 'Q100%, Ca 2+ ' is shown in blue.For visual clarity, in addition to the APL with standard deviation (left), we show the running average (over 5 ns windows) without standard deviation (right).The plots are further split into two showing APL data from simulations with Ca 2+ and Na + neutralizing cations.S11.Schematic illustrations showing the modified Martini 2 bead names and types, and the Lennard-Jones (LJ) interaction matrix compared between the standard and fast Martini 2 ReLPS parameters.(A) ReLPS labelled as per the standard Martini 2 bead names and types.These beads are coloured according to our parametrization of the fast Martini 2 ReLPS.The headgroup beads (green) are named with a prefix 'R' given that they undergo a repulsion of 2 levels compared to standard Martini 2, and all other beads are named 'M' given that they had to be modified considering their potential interactions with the 'R' beads.(B) Comparison of LJ interaction matrices between the standard and fast Martini 2 ReLPS beads.LJ levels 0 to 9 are coloured on a blue-white-red (low-high) scale.The epsilon values of LJ potentials that are altered by 'n' times in the fast ReLPS model with reference to the standard model are also shown in a matrix.S12.Evaluating the influence of the fast Martini 2 ReLPS parameters on its MSD, NNI, and partial densities.(A) Comparing the MSDs and (B) NNI of (1) fast Martini 2 ReLPS (magenta), (2) ReLPS with only their charges modified i.e. 'Q50% Na + ' (orange), (3) ReLPS with only their LJ modified i.e. 'Head LJ+2' (green), and (4) the reference system with standard Martini 2 parameters (blue).(C) Partial density profiles of various components of the system containing fast Martini 2 ReLPS (top half of the plot) compared to those of the standard Martini 2 (unmodified) system (bottom half of the plot).The density profiles are such that the hydrophobic core of the membrane is located at the i.e.Z distance = 0 nm (grey dotted line), and the periplasmic and extracellular regions are located on the left and right hand side of the centre respectively.S13.Comparing APL, membrane thickness, and MSDs of Martini 2 systems containing OMPs surrounded by 100% fast ReLPS, 100% standard ReLPS, and 50% fast/standard ReLPS in the outer leaflet.(A) Averaged APL from both membrane leaflets and (B) membrane thickness versus time of systems containing 100% fast ReLPS (magenta), 100% standard ReLPS (cyan), and 50% fast/standard ReLPS (light red) in the outer leaflet of the outer membrane model.Standard deviations for each data line is shown in light grey.(C) MSDs versus time compared between 100% fast and 100% standard Martini 2 (left).In the 50% fast/standard ReLPS-containing systems, two initial setups were obtained i.e. a non-uniform and a uniform distribution of Martini 2 fast/standard ReLPS.From each of these simulations, MSDs were calculated separately for the fast ReLPS molecules and for the standard ReLPS molecules.A comparison of MSDs between the fast and standard ReLPS molecules is made from the simulations starting with the non-uniform distribution (center) and the uniform distribution (right).(D) Snapshots of the initial setup (left) and final frame (right) of simulations starting with a uniform distribution (top) and non-uniform distribution (bottom) of fast and standard ReLPS molecules.S14.Internal flexibility and lateral mobility of ReLPS in a 100% ReLPS symmetric bilayer.(A,B,C) RMSD histograms for 5 individual ReLPS molecules.The full 10 s trajectory was analyzed for the (A) standard Martini 2 and (B) fast Martini 2 models.RMSD was calculated on the whole molecule after fitting in the central region (GlcN and phosphates).The atomistic simulation was mapped into Martini 2 where 5 randomly chosen individual ReLPS molecules were analyzed as explained and shown in (C).(D) Lateral MSD of ReLPS in replicas of standard (M2 r1-3) and fast Martini 2 (M2 f1-3) ReLPS computed independently for 10 non-overlapping 1 µs-long trajectory segments.S17.DOPC and PMB1 insertion tests conducted with standard and fast Martini 2 ReLPS.(A) Distance along the Z dimension (perpendicular to the membrane) versus time calculated between the phosphate bead(s) of the (single) DOPC molecule and of the ReLPS molecules in the outer leaflet of the membrane.(B) Average number of contacts of PMB1 beads with standard Martini 2 ReLPS beads are mapped onto their respective coarse-grained structures on a blue-white-red (low-high) scale.(C) The same is performed for contacts between PMB1 and fast Martini 2 ReLPS.For both (A) and (B), their respective contact probabilities are plotted on a matrix and shown below using the same colour scale.For visual clarity, labels are alternated between opposite axes.S18.Steered MD simulations of DPC.(A) Snapshots taken from steered MD simulations of DPC in Martini 3. (B) Data showing the change in pulling force versus time while DPC was steered from water into the membrane using Martini 2 (left) and Martini 3 ReLPS (right) models.In each plot, the maximum force experienced by DPC is labelled and corresponds to the yellow region in each plot.S19.Densities of ReLPS relative to OmpF in an OM model with 2% ReLPS in the outer leaflet.(A) Radial distribution functions and (B) 2D density maps of the Martini 2 and 3 ReLPS models with reference to an OmpF trimer.S21.Dissociation events between Martini 2 ReLPS models and OmpF in an OM model with 2% ReLPS in the outer leaflet.Number of OmpF-ReLPS contacts versus time shown from three individual simulations conducted with (A) standard and (B) fast Martini 2 ReLPS.Some dissociation events are indicated by arrows below the respective plots with their colours matching with that of the ReLPS molecule number.S22.Dissociation events between our Martini 3 ReLPS model and OmpF in an OM model with 2% ReLPS in the outer leaflet.Number of OmpF-ReLPS contacts versus time shown from three individual simulations conducted with our Martini 3 ReLPS model.Some dissociation events are indicated by arrows below the respective plots with their colours matching with that of the ReLPS molecule number.

Table S3 .
Choosing an LJ-modified model of Martini 2 ReLPS such that the difference between the APL of the standard Martini 2 (Head LJ+0) and of the modified Martini 2 model (Head LJ+n) is < 5% (n>0).

Table S4 .
Details on box dimensions per system and GROMACS versions to perform production simulations.Box dimension units are in (X x Y x Z) nm.Note that the box dimensions are shown for 3 replicas at the end of their production simulations, including those systems which were simulated for more than 3 replicas.In the third column, all box dimensions that were setup before energy minimization are mentioned, except the cells where "eq0" is mentioned in which case box dimensions were obtained after a short equilibration.Abbreviations used in the table are: C36m: CHARMM36m, M2: Martini 2 standard, M2 fast: fast Martini 2, M3: Martini 3, OM: outer membrane.